pomp
object,
encoding a partially-observed Markov process model together with a uni- or multi-variate time series.
As such, it is central to all the package's functionality.
One implements the model by specifying some or all of its basic components.
These include:
The basic structure and its rationale are described in the Journal of Statistical Software paper, an updated version of which is to be found on the package website.
pomp(data, times, t0, ..., rprocess, dprocess, rmeasure, dmeasure, measurement.model, skeleton, skeleton.type, skelmap.delta.t, initializer, rprior, dprior, params, covar, tcovar, obsnames, statenames, paramnames, covarnames, zeronames, PACKAGE, fromEstimationScale, toEstimationScale, globals, cdir, cfile, shlib.args)
data
should be given as a data-frame and times
must indicate the column of observation times by name or index.
times
must be numeric and strictly increasing.
Internally, data
will be internally coerced to an array with storage-mode double
. In addition, a pomp
object can be supplied in the data
argument.
In this case, the call to pomp
will add element to, or replace elements of, the supplied pomp
object.
t0 <= times[1]<="" code="">.
=>
Note: it is not typically necessary (or even feasible) to define dprocess
.
In fact, no current pomp inference algorithm makes use of dprocess
.
This functionality is provided only to support future algorithm development.
params
and an initial time, t0
, the initializer determines the state vector at time t0
.
See below under The State-Process Initializer for details.
double
.
covar
is the table of covariates (one column per variable);
tcovar
the name or the index of the time variable. If a covariate table is supplied, then the value of each of the covariates is interpolated as needed.
The resulting interpolated values are made available to the appropriate basic components.
Note that covar
will be coerced internally to storage mode double
.
See below under Covariates for more details.
obsnames
or covarnames
, as these will by default be read from data
and covars
, respectively.
toEstimationScale
and fromEstimationScale
are transformations from the model scale to the estimation scale, and vice versa, respectively.
See below under Parameter Transformations for more details.
pomp
uses C snippets.
If no C snippets are used, globals
has no effect.
cdir
specifies the name of the directory within which C snippet code will be compiled.
By default, this is in a temporary directory specific to the running instance of R.
cfile
gives the name of the file (in directory cdir
) into which C snippet codes will be written.
By default, a random filename is used.
The shlib.args
can be used to pass command-line arguments to the R CMD SHLIB
call that will compile the C snippets.
pomp
will be made available to each of the basic components.
To prevent errors due to misspellings, a warning is issued if any such arguments are detected.
pomp
constructor function returns an object, call it P
, of class pomp
.
P
contains, in addition to the data, any elements of the model that have been specified as arguments to the pomp
constructor function.
One can add or modify elements of P
by means of further calls to pomp
, using P
as the first argument in such calls.
pomp
constructor causes them to be written to a C file stored in the R session's temporary directory, which is then compiled (via R CMD SHLIB
) into a dynamically loadable shared object file.
This is then loaded as needed. Note to Windows and Mac users:
By default, your R installation may not support R CMD SHLIB
.
The package website contains installation instructions that explain how to enable this powerful feature of R.R CMD SHLIB
.
If the resulting file does not compile, an error message will be generated.
No attempt is made by pomp to interpret this message.
Typically, compilation errors are due to either invalid C syntax or undeclared variables.
statenames
or paramnames
arguments to pomp
, respectively.
Compiler errors that complain about undeclared state variables or parameters are usually due to failure to declare these in statenames
or paramnames
, as appropriate.
rmeasure
and dmeasure
) by those names.
pomp
object contains a table of covariates (see above), then the variables in the covariate table will be available, by their names, in the context within which the C snippet is executed.
file.show(system.file("include/pomp.h",package="pomp"))to view this header file.
globals
argument of pomp
will be included at the head of the generated C file.
This can be used to declare global variables, define useful functions, and include arbitrary header files.
rprocess
and/or dprocess
is facilitated by pomp's so-called plug-ins, which allow one to easily specify the most common kinds of process model. Discrete-time processes
If the state process evolves in discrete time, specify rprocess
using the discrete.time.sim
plug-in.
Specifically, provide rprocess = discrete.time.sim(step.fun = f, delta.t)to
pomp
, where f
is a C snippet or R function that takes simulates one step of the state process.
The former is the preferred option, due to its much greater computational efficiency.
The goal of such a C snippet is to replace the state variables with their new random values at the end of the time interval.
Accordingly, each state variable should be over-written with its new value.
In addition to the states, parameters, covariates (if any), and observables, the variables t
and dt
, containing respectively the time at the beginning of the step and the step's duration, will be defined in the context in which the C snippet is executed.
See below under General rules for C snippet writing for more details.
Examples are to be found in the tutorials on the package website. If f
is given as an R function, it should have prototype f(x, t, params, delta.t, ...)When
f
is called, x
will be a named numeric vector containing the value of the state process at time t
,
params
will be a named numeric vector containing parameters,
and delta.t
will be the time-step.
It should return a named vector of the same length, and with the same set of names, as x
, representing a draw from the distribution of the state process at time t+delta.t
, conditional on its having value x
at time t
. Continuous-time processes
If the state process evolves in continuous time, but you can use an Euler approximation, specify rprocess
using the euler.sim
plug-in.
Furnish rprocess = euler.sim(step.fun = f, delta.t)to
pomp
in this case.
As before, f
can be provided either as a C snippet or as an R function, the former resulting in much quicker computations.
The form of f
will be the same as above (in the discrete-time case). If you have a procedure that allows you, given the value of the state process at any time, to simulate it at an arbitrary time in the future, use the onestep.sim
plug-in.
To do so, furnish rprocess = onestep.sim(step.fun = f)to
pomp
.
Again, f
can be provided either as a C snippet or as an R function, the former resulting in much quicker computations.
The form of f
should be as above (in the discrete-time or Euler cases). If you desire exact simulation of certain continuous-time Markov chains, an implementation of Gillespie's algorithm (Gillespie 1977) is available, via the gillespie.sim
plug-in.
In this case, furnish rprocess = gillespie.sim(rate.fun = f, v, d)to
pomp
, where f
gives the rates of the elementary events.
Here, f
must be an R function of the form f(j, x, t, params, ...)When
f
is called,
the integer j
will be the number of the elementary event (corresponding to the columns of matrices v
and d
, see below),
x
will be a named numeric vector containing the value of the state process at time t
and
params
is a named numeric vector containing parameters.
f
should return a single numerical value, representing the rate of that elementary event at that point in state space and time. Matrices v
and d
specify the continuous-time Markov process in terms of its elementary events.
Each should have dimensions nvar
x nevent
, where nvar
is the number of state variables and nevent
is the number of elementary events.
v
describes the changes that occur in each elementary event:
it will usually comprise the values 1, -1, and 0 according to whether a state variable is incremented, decremented, or unchanged in an elementary event.
d
is a binary matrix that describes the dependencies of elementary event rates on state variables:
d[i,j]
will have value 1 if event rate j
must be updated as a result of a change in state variable i
and 0 otherwise. A faster, but approximate, version of the Gillepie algorithm, the so-called K-leap method of Cai and Xu (2007), is implemented in the kleap.sim
plug-in.
To use it, supply rprocess = kleap.sim(rate.fun, e, v, d)to
pomp
, where rate.fun
, v
, and d
are as above, and e
gives relative error tolerances for each of the state variables.
e
should have length equal to the number of state variables.
e[i]
corresponds to row i
of the v
and d
matrices and we must have $0 <= e[i]="" <="1$." the="" leap="" size,="" $k$,="" is="" chosen="" so="" that="" $k="" size="" of="" time="" step="" simulator="" plug-ins="" discrete.time.sim, euler.sim
, and onestep.sim
all work by taking discrete time steps.
They differ as to how this is done.
Specifically,
-
onestep.sim
takes a single step to go from any given time t1
to any later time t2
(t1 < t2
).
Thus, this plug-in is designed for use in situations where a closed-form solution to the process exists.
- To go from
t1
to t2
, euler.sim
takes n
steps of equal size, where
n = ceiling((t2-t1)/delta.t).
-
discrete.time.sim
assumes that the process evolves in discrete time, where the interval between successive times is delta.t
.
Thus, to go from t1
to t2
, discrete.time.sim
takes n
steps of size exactly delta.t
, where
n = floor((t2-t1)/delta.t).
Specifying dprocess
If you have a procedure that allows you to compute the probability density of an arbitrary transition from state $x1$ at time $t1$ to state $x2$ at time $t2$, assuming that the state remains unchanged between $t1$ and $t2$, then you can use the onestep.dens
plug-in.
This is accomplished by furnishing
dprocess = onestep.dens(dens.fun = f)
to pomp
, where f
is an R function with prototype
f(x1, x2, t1, t2, params, ...)
When f
is called,
x1
and x2
will be named numeric vectors containing the values of the state process at times t1
and t2
, respectively,
and params
will be a named numeric vector containing parameters.
f
should return the log likelihood of a transition from x1
at time t1
to x2
at time t2
, assuming that no intervening transitions have occurred. To see examples, consult the tutorials on the package website. =>
rprocess
and dprocess
arguments, or via the measurement.model
argument.
If measurement.model
is given it overrides any specification via the rmeasure
or dmeasure
arguments, with a warning. The best way to specify the measurement model is by giving C snippets for rmeasure
and dmeasure
.
In writing an rmeasure
C snippet, bear in mind that:
t
, containing the time of the observation, will be defined in the context in which the snippet is executed.
rmeasure
using an R function.
In this case, specify the measurement model simulator by furnishing rmeasure = fto
pomp
, where f
is an R function with prototype f(x, t, params, \dots)It can also take any additional arguments if these are passed along with it in the call to
pomp
.
When f
is called,
x
will be a named numeric vector of length nvar
, the number of state variables.
t
will be a scalar quantity, the time at which the measurement is made.
params
will be a named numeric vector of length npar
, the number of parameters.
f
must return a named numeric vector of length nobs
, the number of observable variables. In writing a dmeasure
C snippet, observe that:
t
, containing the time of the observation,
and the Boolean variable give_log
will be defined in the context in which the snippet is executed.
lik
variable to the likelihood of the data given the state.
Alternatively, if give_log == 1
, lik
should be set to the log likelihood.
dmeasure
is to be provided instead as an R function, this is accomplished by supplying dmeasure = fto
pomp
, where f
has prototype f(y, x, t, params, log, \dots)Again, it can take additional arguments that are passed with it in the call to
pomp
.
When f
is called,
y
will be a named numeric vector of length nobs
containing values of the observed variables;
x
will be a named numeric vector of length nvar
containing state variables;
params
will be a named numeric vector of length npar
containing parameters;
t
will be a scalar, the corresponding observation time.
f
must return a single numeric value, the probability density of y
given x
at time t
.
If log == TRUE
, then f
should return instead the log of the probability density.
Note: it is a common error to fail to account for both log = TRUE
and log = FALSE
when writing the dmeasure
C snippet or function. One can also specify both the rmeasure
and dmeasure
components at once via the measurement.model
argument.
It should be a formula or list of nobs
formulae.
These are parsed internally to generate rmeasure
and dmeasure
functions.
Note: this is a convenience function, primarily designed to facilitate model exploration;
it will typically be possible (and as a practical matter necessary) to accelerate measurement model computations by writing dmeasure
and/or rmeasure
using C snippets.trajectory
and traj.match
. If the state process is a discrete-time stochastic process, then the skeleton is a discrete-time map.
To specify it, provide skeleton = map(f, delta.t)to
pomp
, where f
implements the map and delta.t
is the size of the timestep covered at one map iteration. If the state process is a continuous-time stochastic process, then the skeleton is a vectorfield (i.e., a system of ordinary differential equations).
To specify it, supply skeleton = vectorfield(f)to
pomp
, where f
implements the vectorfield, i.e., the right-hand-size of the differential equations. In either case, f
can be furnished either as a C snippet (the preferred choice), or an R function.
In writing a skeleton
C snippet, be aware that:
x
is named Dx
and is the new value of x
after one iteration of the map.
x
is named Dx
and is the value of $dx/dt$.
t
, will be defined in the context within which the snippet is executed.
f
is an R function, it must be of prototype f(x, t, params, \dots)where, as usual,
x
is a numeric vector (length nvar
) containing the coordinates of a point in state space at which evaluation of the skeleton is desired.
t
is a scalar value giving the time at which evaluation of the skeleton is desired.
params
is a numeric vector (length npar
) holding the parameters.
f
may take additional arguments, provided these are passed along with it in the call to pomp
.
The function f
must return a numeric vector of the same length as x
, which contains the value of the map or vectorfield at the required point and time.t0
).
By default, pomp
assumes that this initial distribution is concentrated on a single point.
In particular, any parameters in params
, the names of which end in .0
, are assumed to be initial values of states.
When the state process is initialized, these are simply copied over as initial conditions.
The names of the resulting state variables are obtained by dropping the .0
suffix. One can override this default behavior by furnishing a value for the initializer
argument of pomp
.
As usual, this can be provided either as a C snippet or as an R function.
In the former case, bear in mind that:
t
, containing the zero-time, will be defined in the context in which the snippet is executed.
statenames
argument plays a particularly important role when the initializer is specified using a C snippet.
In particular, every state variable must be named in statenames
.
Failure to follow this rule will result in undefined behavior.
initializer = fto
pomp
, where f
is a function with prototype f(params, t0, \dots)When
f
is called,
params
will be a named numeric vector of parameters.
t0
will be the time at which initial conditions are desired.
f
may take additional arguments, provided these are passed along with it in the call to pomp
.
f
must return a named numeric vector of initial states.
It is of course important that the names of the states match the expectations of the other basic components. Note that the state-process initializer can be either deterministic (the default) or stochastic.
In the latter case, it samples from the distribution of the state process at the zero-time, t0
.rprior
and/or dprior
arguments to pomp
.
As with the other basic model components, it is preferable to specify these using C snippets.
In writing a C snippet for the prior sampler (rprior
), keep in mind that:
dprior
), observe that:
give_log
will be defined.
lik
.
When give_log == 1
, lik
should contain the log of the prior probability density.
rprior
must be a function of prototype f(params, \dots)that makes a draw from the prior distribution given
params
and returns a named vector of the same length and with the same set of names, as params
.
The dprior
function must be of prototype f(params, log = FALSE, \dots).Its role is to evaluate the prior probability density (or log density if
log == TRUE
) and return that single scalar value.pomp
object contains covariates (specified via the covar
argument; see above), then interpolated values of the covariates will be available to each of the model components whenever it is called.
In particular, variables with names as they appear in the covar
data frame will be available to any C snippet.
When a basic component is defined using an R function, that function will be called with an extra argument, covars
, which will be a named numeric vector containing the interpolated values from the covariate table. An exception to this rule is the prior (rprior
and dprior
): covariate-dependent priors are not allowed.
Nor are parameter transformations permitted to depend upon covariates.pomp
object via the toEstimationScale
and fromEstimationScale
arguments.
As with the basic model components, these should ordinarily be specified using C snippets.
When doing so, note that:
toEstimationScale
argument whilst the transformation mapping a parameter vector from the alternative scale to that on which the model is defined is specified with the fromEstimationScale
argument.
p
should be assigned to variable Tp
.
params
and ...
.
In this case, toEstimationScale
should transform parameters from the scale that the basic components use internally to the scale used in estimation.
fromEstimationScale
should be the inverse of toEstimationScale
. Note that it is the user's responsibility to make sure that the transformations are mutually inverse.
If obj
is the constructed pomp
object, and coef(obj)
is non-empty, a simple check of this property is x <- coef(obj, transform = TRUE) obj1 <- obj coef(obj1, transform = TRUE) <- x identical(coef(obj), coef(obj1)) identical(coef(obj1, transform=TRUE), x)By default, both functions are the identity transformation.
pomp
's zeronames
argument will be set to zero immediately following each observation.
See euler.sir
and the tutorials on the package website for examples.pomp
with one or more C snippet arguments.
You can set cdir
and cfile
to control where this code is written.
Alternatively, set options(verbose=TRUE)
before calling pomp
.
This will cause a message giving the name of the generated C file (in the session temporary directory) to be printed.pomp
, but complete error checking for arbitrary models is impossible.
If the user-specified functions do not conform to the above specifications, then the results may be invalid.
In particular, if both rmeasure
and dmeasure
are specified, the user should verify that these two functions correspond to the same probability distribution.
If skeleton
is specified, the user is responsible for verifying that it corresponds to a deterministic skeleton of the model.D. T. Gillespie (1977) Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry 81:2340--2361. X. Cai and Z. Xu (2007) K-leap method for accelerating stochastic simulation of coupled chemical reactions. Journal of Chemical Physics 126:074102.
## pomp encoding a stochastic Ricker model with a covariate:
pomp(data = data.frame(t = 1:100, y = NA),
times = "t", t0 = 0,
covar = data.frame(t=0:100,K=seq(from=50,to=200,length=101)),
tcovar = "t",
rprocess = discrete.time.sim(Csnippet("double e = rnorm(0,sigma);
n = r*n*exp(1-n/K+e);"), delta.t = 1),
skeleton = map(Csnippet("Dn = r*n*exp(1-n/K);"), delta.t = 1),
rmeasure = Csnippet("y = rpois(n);"),
dmeasure = Csnippet("lik = dpois(y,n,give_log);"),
rprior = Csnippet("r = rgamma(1,1); sigma = rgamma(1,1);"),
dprior = Csnippet("lik = dgamma(r,1,1,1) + dgamma(sigma,1,1,1);
if (!give_log) lik = exp(lik);"),
initializer = Csnippet("n = n_0;"),
toEstimationScale = Csnippet("Tr = log(r); Tsigma = log(sigma);"),
fromEstimationScale = Csnippet("Tr = exp(r); Tsigma = exp(sigma);"),
paramnames = c("n_0", "r", "sigma"),
statenames = "n") -> rick
## fill it with simulated data:
rick <- simulate(rick, params = c(r=17, sigma = 0.1, n_0 = 50))
plot(rick)
## Not run:
# pompExample()
# demos(package="pomp")
# ## End(Not run)
Run the code above in your browser using DataLab